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The diffusive characteristics of two upwind schemes, multi-dimensional fluctuation 
splitting and dimensionally-split finite volume, are compared for scalar advection-diffusion 
problems. Algorithms for the two schemes are developed for node-based data represen- 
tation on median-dual meshes associated with unstructured triangulations in two spatial 
dimensions. Four model equations are considered: linear advection, non-linear advection, 
diffusion, and advection-diffusion. Modular coding is employed to isolate the effects of the 
two approaches for upwind flux evaluation, allowing for head-to-head accuracy and effi- 
ciency comparisons. Both the stability of compressive limiters and the amount of artificial 
diffusion generated by the schemes is found to be grid-orientation dependent, with the 
fluctuation splitting scheme producing less artificial diffusion than the dimensionally-split 
finite volume scheme. Convergence rates are compared for the combined advection-dif- 
fusion problem, with a speedup of 2—3 seen for fluctuation splitting versus finite volume 
when solved on the same mesh. However, accurate solutions to problems with small 
diffusion coefficients can be achieved on coarser meshes using fluctuation splitting rather 
than finite volume, so that when comparing convergence rates to reach a given accuracy, 
fluctuation splitting shows a 20—25 speedup over finite volume. 


Nomenclature 

A, B Flux Jacobians 

c Nodal update coefficients 

F Convective flux vector, F = F(x, y,u ) 
h Mesh spacing 

i, j Cartesian unit vectors 

J -1 Inverse Jacobian of the coordinate transfor- 
mation 

£ Edge length 

h Outward unit normal to control cell 

n Length-scaled inward normal 

Q Fluctuation ratio 

r Distance vector from node to face 

Si Median-dual area about node i 
t Time 

u Dependent variable 

x,y Cartesian coordinates 

a, 3 Curvilinear advection speeds 

, Boundary of control cell 

7 Limiter bound 

d Finite element linear shape function 

A Advection velocity vector 

v Diffusion coefficient 
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£, g Curvilinear coordinates 

r Timestep 

$ Artificial dissipation — finite volume 

< p Element advective fluctuation 

( ft , Fluctuation components 

<p ri Limited fluctuations 

(jfi , <A Artificial dissipation — fluctuation splitting 
<j> v Viscous fluctuation 

tb Limiter function 

D Area of control cell 

V Gradient operator 

Introduction 

U PWIND discretizations for advection equations 
typically introduce artificial numerical dissipa- 
tion into the solution. When combined advection- 
diffusion problems are considered, this dissipation in- 
troduced in the discretization of the advection terms 
should be less than the true physical diffusion. To 
this end the diffusive characteristics of upwind schemes 
are investigated and their performance in resolving 
solutions to advection-diffusion problems with small 
diffusion coefficients is analyzed. 

Two node-based, median-dual methods for modeling 
convective fluxes are considered. The first is a tradi- 
tional dimensionally- split finite volume (FV) scheme. 1 
Dimensionally-split schemes are known to introduce 
excess dissipation when discontinuities are not aligned 
with the mesh. 2 

The second method is the NNL 3 fluctuation splitting 
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(FS) scheme. FS has a more-compact stencil than FV 
for second-order formulations and exhibits “zero cross- 
diffusion”! in a grid-aligned condition. Both of these 
attributes should lead to less introduced dissipation 
as compared with FV. The sensitivity of FS to grid 
orientation and resulting production of cross-diffusion 
is investigated in the present report. 

The use of compressive limiter functions is inves- 
tigated with both algorithms. Also, local timesteps 
based on positivity arguments are tested, and in some 
cases the optimum timestep for convergence is found 
to be less than the largest stable timestep. This be- 
havior is attributed to residual oscillations forced by 
“ringing” of the limiter functions. 

Formulation of FS schemes for diffusion problems 
is a recent research area. 4,5 The present study seeks 
to quantify the relative merits of using a low-diffusion 
advection operator to resolve advection-diffusion prob- 
lems with small diffusion coefficients. Lessons learned 
on these problems will guide the development of com- 
puter codes for solving compressible viscous fluid dy- 
namic problems. 

Governing Equations 

The non-linear advection/diffusion equation, 

u t + V-F = V-(i/Vu) (1) 

is cast as a hyperbolic conservation law, to which 
steady-state solutions are sought. 

Finite Volume 

In FV form, using the divergence theorem Eqn. 1 
becomes, 


/ 

J n 


ut clfl 



— i /Vu) ■ n d. 


( 2 ) 


where SI is the median dual about node i and , is 
the boundary of fl. Using mass lumping to the nodes, 
similar to an explicit finite element treatment, 6 the 
temporal evolution is evaluated on a time-invariant 
mesh as, 


I Ut dn [u'r - u';) (3) 

The discretization of the convective flux, F, is per- 
formed using Barth’s implementation 1 of the upwind, 
locally one-dimensional, approximate Riemann prob- 
lem of Roe. 7 


® F ■ n d, -A ^ 

2 {^in “1“ Font) ’ ^ ^ 

’ faces 

L J 


A, (4) 


t “Zero cross diffusion” refers to the practice of adding artifi- 
cial diffusion terms in the streamwise direction only, as opposed 
to adding artificial dissipation in both the streamwise and cross- 
stream directions. 


where the artificial dissipation provides the upwinding, 

T = )fAn x . + Bn,j\(u 0Ut - u in ) (5) 

with n = n x i + n v j. Out and in refer to states on the 
outside and inside of 0 at the face. A and B are the 
flux Jacobians, 


apd) 3p(2) 

' “ ~ihT : B - ~^T 


( 6 ) 


and ( A,B ) represent their conservative linearizations 
at the cell face. 7 

Reconstruction from the nodal unknowns to the cell 
faces, as, 


Uface = Ui + V’Vu • T (7) 

is performed with Barth’s limited, unweighted least- 
squares procedure to provide second-order spatial ac- 
curacy in smoothly- varying regions of the solution. 

Two methods for evaluating the diffusion term are 
incorporated into FV. The more compact of the two, 
the finite element discretization, is discussed in the 
following section. The less-compact diffusion formula 
is obtained by discretizing the last term of Eqn. 2, in 
a manner similar to Eqn. 4, 

j ( v V u ■ h d, — > 

'U'in + Vu 0Ut )-n A, (8) 

’ faces 

The diffusion coefficient is averaged over the length of 
the face. The gradients from Eqn. 7 are not limited 
before averaging at the control-volume faces in Eqn. 8, 
as suggested by Anderson and Bonhaus. 8 


Fluctuation Splitting 

The NNL FS scheme is presented as a slight re- 
interpretation of the work of Sidilkover and Roe. 2 3 The 
current interpretation is as a volume integral over tri- 
angular elements, without recourse to the divergence 
theorem. The discretized equations, however, are iden- 
tical. 

Integrating Eqn. 1 over an element, where Q is now 
the area of the triangular element, 

[ u t dn = - [ V • F dn + [ V • (i/Vm) dn (9) 
J Q J Q J Q 

For linear variation of the dependent variable over the 
element, the temporal evolution is, 


/ 

Jn 


Ut dn = f lu t = + u '2 , + u 3t ) 


( 10 ) 


where u\, u 3 , and u 3 correspond to the three nodes 
defining element fl. 

Defining local curvilinear coordinates, £ and g, par- 
allel to sides 12 and 23, respectively (Fig. 1), the 
divergence of the convective flux can be written, 

V • F = F^ + Ef = (h 2 • Ft - hi • Fri) (11) 
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3 


with, 



Fig. 1 Fluctuation splitting element nomencla- 
ture. 


Defining the scaled inward normal, n = —hn, where 
h is a mesh edge length, the divergence (Eqn. 11) be- 
comes, 


— ( 
20 V 


5 + ^23 n l ' (12) 


V • F = — -h 12 n 2 ■ F t 


If F is linear or quadratic in u, then for a linear 
variation of u over the element, 


/ 

Jn 


V • F dfl = aA 2 iu -1- / 3A 22 u (13) 


where the difference operator is defined A 2 iu = u 2 —u\ 
and the advection speeds are, 

a = — ^( n 2,d + n2 s B), /3 = ^(n^I + nqB) (14) 

A and B are now the conservative linearizations over 
the triangular element. 9 

The advective fluctuation can be defined. 


(/> = -[ v • f (in 

Jn 


(15) 


(16) 


The fluctuation can be split, 
where, 

= oAjim, o'' = -3A S2 u (17) 

Following Sidilkover 10 the fluctuation is limited to 
achieve a second-order scheme, 

d> r = ^ + <p n ip{Q) = 0^ ^1 - yy j (18) 


Q = - 


( i> v 


(20) 


Upwinding is achieved through the introduction of 
the artificial dissipation terms, 

<jj ^ = sign (a)4>^ , (f> v = sign(/3)<)>’ , (21) 


Combining Eqn. 10 with a distribution scheme for 
Eqn. 15 and summing over all elements, the contribu- 
tions to nodal time derivatives can be written in the 
form, 


.SV«i, 

<r~ 


-<!>*)+ COE 

SoU-2 t 

F- 




1 


S 2 U3 t 

1- 

2^ 

+ <t> ri )+COE 


( 22 ) 


where COE stands for contributions from other ele- 
ments containing these nodes. 

A finite element treatment, similar to Tomaich, 3 4 is 
employed to obtain the diffusive fluctuation, 


4>v = / V • (i/Vm) 

Jn 


dfl 


(23) 


Assuming piecewise-linear data and an element-aver- 
aged diffusion coefficient leads to a diffusive fluctuation 
of zero for the triangular element. Introducing the 
linear nodal shape functions such that = 1> 

the elemental diffusive fluctuation can be expressed 
< i>v yi ! , (pvi = o, where 


4> Vi = I d, : V • (vVu) 

Jn 


dfl 


(24) 


Integrating by parts, 

4> Vi = j dji'Y ii • fid , — j vS!u ■ Vd,: dfl (25) 

The boundary integral in Eqn. 25 will cancel on sum- 
ming contributions for interior nodes. The remaining 
volume integral can be evaluated analytically, 

3 

<p Vi = TjV,/ • 11 ;. | = — ^2 "J"J ' ■ ' ■ ( 26 ) 

J = 1 

Distributing this diffusive fluctuation to the nodes and 
keeping only the larger of the physical or artificial dis- 
sipation leads to the update formula, 

Siu u <- y + max|-y,^ij +COE 

Sam, m A±fL + max ( <£^£l,, h 3 + coe 
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Boundary Conditions 

Explicit Dirichlet inflow boundary conditions are 
employed. Advective outflow boundaries are treated 
for free convection through the boundary nodes, al- 
lowing boundary nodes to be handled in the same 
manner as interior nodes. For the diffusion terms a 
Neumann outflow boundary is applied with zero gra- 
dient, achieved by setting the boundary integral in 
Eqn. 25 to zero. 


Limiter Functions 

The limiter functions employed by both schemes are 
cast as a limit on the ratio of two values. The limiting 
proposed by Barth 1 is interpreted as a non-symmetric 
limiter 11 consistent with one-dimensional limiting in 
the maximum gradient direction. The limiting on the 
reconstruction (Eqn. 7) is performed on the ratio, 


ip 


il> 


yVninfmax 

2 Vw, • r 


(28) 


where u min / max is the minimum (resp. maximum) of 
Uj and all distance-one neighbors. This same inter- 
pretation has been put forth by Bruner and Wal- 
ters , 12 though they identified Barth’s limiting with 
the symmetric Superbee limiter, rather than the non- 
symmetric Cliakravarthy and Osher 11 limiter. 

Minmod, van Albada , 13 Superbee, and 7 14 sym- 
metric limiters are utilized for FV (Eqn. 7) and FS 
(Eqns. 18 and 19). The 7 limiter, of which the Min- 
mod (7 = 1) and Superbee (7 = 2) are special cases, 
is, 


tp'y(Q) = max[0, inin( 7 < 3 ,l), min(Q, 7 )] (29) 


Timestep 

Both schemes are formulated either as Gauss-Seidel 
time-relaxation or forward Euler time-evolution algo- 
rithms. 

The nodal updates for the discrete system can be 
formed as a sum of contributions from all nodes. 

u) +r = r i"j = r ;"i + r '"./ ( 30 ) 

i 

For positivity 15 each of the coefficients in Eqn. 30 must 
be non-negative. 

Advective Timestep restriction 

In the FV context the nodal update (Eqn. 30) can 
be rearranged into the form of Eqn. 3, 

y(u i +r - u\) = 77 (r '- - I)"/ c i"j ( 31 ) 

iAi 

For the upwind, edge-based algorithm considered here, 
each ^-Cj will be related to a positive-definite coeffi- 
cient equal to zero for outflowing faces and related to 


the wavespeed for inflowing faces, yielding the restric- 
tion t > 0 on the timestep. The remaining term can 
be expressed, 

~( c i ~ 1) = - E c * ( 32 ) 

k about i 

where the c & coefficients are also positive-definite, ei- 
ther zero for inflowing faces or related to the wavespeed 
for outflowing faces. Rearranging and imposing the 
positivity constraint, c, > 0 , yields the timestep re- 
striction, 

1 - 7 T y> = c,>0 (33) 

k about l 


r < - - (34) 

Z^k about i ^k 

For FS, the nodal updates are assembled from 
Eqn. 22 as, 

^(u* +r - u\) = Cj ( Uj - m) (35) 

In this case the cj coefficients are formed as contri- 
butions from the fluctuations in the triangles to both 
the left and the right of mesh edge rj. The positivity 
restriction on r is found to have a similar form as for 
finite volume (Eqn. 34), 

r < (36) 

c -> 

Local time-stepping based on positivity is shown 
to yield stable, yet non-converging, solutions in some 
second-order cases (see Results section). Robust con- 
vergence is obtained by using the first-order c’s in 
Eqns. 34 and 36, even for second-order-accurate spa- 
tial discretizations. 


Diffusive Timestep Restriction 

Unfortunately, the finite element formulation for the 
diffusive terms (Eqn. 26) cannot be guaranteed to pre- 
serve local positivity on obtuse triangles (see Barth 1 ). 
Considering only the contributions from the current 
node, the coefficient for the diffusion term can be writ- 
ten, 



The appropriate edge length is the side of the element 
opposite the current node. The resulting timestep re- 
striction is, 


r < 


vt 2 

40 


(38) 


In a similar manner the timestep restriction from 
Eqn. 8 is, 


r < 


Si 


3zzA, 2 
2-^T 4Q 


( 39 ) 
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Fig. 2 First-order fluctuation splitting, uniform 
advection. 

Results 

Linear Advection 

The linear advection equation is obtained from 
Eqn. 1 by setting v = 0 and F = Xu, yielding, 

u t + V • (Xu) = 0 (40) 

A divergence-less advection velocity is considered, such 
that V • A = 0. Equation 40 can then be written, 

ut + A • Vu = 0 (41) 

Uniform Advection 

Uniform advection of the Heavyside function at 
—45 degrees, A = (1,-1), on a cut-cartesian mesh is 
shown for first-order FS, second-order FS, and second- 
order FV in Figs. 2-4, respectively. The mesh is 
shown as the dashed background, and equally-spaced 
contours vary on [0,1], the minimum and maximum 
solution values. The spread of the contour lines with 
spatial evolution is indicative of the amount of dissipa- 
tion introduced into the solution by the discretization 
of the convective terms. 

Second-order FS is seen to be greatly superior to 
first-order, as expected, reproducing the exact solu- 
tion in this case with no introduced dissipation. Also, 
FS represents a significant reduction in numerical dif- 
fusion versus the corresponding FV scheme, with both 
results employing the Minmod limiter. 

However, the “zero cross-diffusion” results of Fig. 3 
with FS are misleading. In Fig. 5 the advection veloc- 
ity has been rotated counter clockwise by 90 degrees 
on the same grid. Clearly, the artificial dissipation 
introduced by the FS scheme has been increased. 

The corresponding F\ T solution is shown in Fig. 6. 
While the change in contour spreading for the FV 
scheme between Figs. 4 and 6 is less dramatic than 
the change in spreading for the FS scheme in Figs. 3 

5 



Fig. 3 Second-order fluctuation splitting, uniform 
advection. 



Fig. 4 Second-order finite volume, uniform advec- 
tion. 

and 5, the FS results still exhibit less diffusion than 
the FV results, comparing Figs. 5 and 6. 

Employing the compressive Superbee limiter with 
the FS scheme yields the results of Fig. 7. In this case 
the discontinuity is confined to a 2-3 cell stencil, and 
does not grow in space. Applying the Superbee lim- 
iter to FV cannot eliminate all artificial dissipation for 
this case, as is possible with FS. The FV results (not 
shown) corresponding to Fig. 7 spread the disconti- 
nuity over four cells by the outflow boundary, with a 
continually broadening trend. 

However, while it is possible to use the Superbee 
limiter with FS for this case, compressive limiters can 
be unstable on different grid orientations. For exam- 
ple, no degree of compression is stable for the case of 
Fig. 3. This potential for instability is related to global 
positivity, as discussed by Sidilkover and Roe. 3 

The effect of using a general unstructured grid is 
investigated in Figs. 8 and 9. The unstructured grid 
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>.=( 1 . 1 ) 

Contour 
spacing: 
0 . 05 - 0.95 
0.1 increments 


Fig. 5 Second-order fluctuation splitting, uniform Fig. 8 Fluctuation splitting on unstructured mesh, 
advection. 



Fig. 9 Finite volume on unstructured mesh. 


Fig. 6 Second-order finite volume, uniform advec- 
tion. 



*. = 0,1 ) 

Contour 
spacing: 
0 . 05 - 0.95 
0.1 increments 


Fig. 7 Second-order fluctuation splitting with 
compressive limiter. 


in this case was generated using Vcrid. 16,17 The FS 
solution exhibits less dissipation, but is not as smooth 
as the FV solution. While the FS scheme preserves 
contact discontinuities over larger spatial ranges than 
the FV scheme, FS does not appear to degenerate 
gracefully with regard to extreme coarsening of the 
unstructured mesh for this test case. This behavior 
could have negative implications for applications em- 
ploying multigrid convergence acceleration. 

Circular Advection 

Circular advection is achieved by setting 
A = (y, —x). A decaying sine- wave input profile 
is used, 

u(x, 0) = (e' r sin7nc) 2 

Results for the two schemes, using the Minmod 
limiter, are presented on the worse-case cut-cartesian 
mesh in Figs. 10 and 11. Again, the FS results are con- 
siderably less diffusive than the dimensionally-split FV 
solution. 
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* = (y. -x) 

Contour 

spacing: 

0.025 - 0.375 
0.05 increments 



l=(y, -x) 

Contour 

spacing: 

0.025 - 0.375 
0.05 increments 


Fig. 10 Fluctuation splitting, circular advection. Fig. 12 Fluctuation splitting on unstructured 

mesh, circular advection. 




Contour 

spacing: 

0.025 - 0.375 
0.05 increments 


U = (e 2x sin(27tx)) i 


Fig. 11 Finite volume, circular advection. 


The circular-advection problem is also applied on 
an unstructured mesh. The input profile for this case 
consists of both a top-hat function and a decaying sine 
wave, allowing comparisons between the schemes for 
both sharp discontinuities and smooth gradients. The 
input profile is, 


u(x , 0) 


(e 2,c sin(27rx))“ —0.5 < x < 0 

0 —0.6 < x < —0.5 

0.4 -0.8 < x < -0.6 

0 -1 < x < -0.8 


Results for this case are displayed in Fig. 12 for FS 
and Fig. 13 for FV, both using the Minmod limiter. FS 
performs significantly better at preserving the top-hat 
distribution. FS also does a better job of maintaining 
the minimum and maximum values of the sine distri- 
bution, though both schemes do well on the smooth 
gradient portion of the sine wave. 


Fig. 13 Finite volume on unstructured mesh, cir- 
cular advection. 

Non-linear Advection 

The non-linear advection equation is obtained from 
Eqn. 1 by setting F = (^, u) with v = 0. In non- 
conservative form the equation is written, 

Ui + ii u. t + u u = 0 

A coalescing shock problem is considered, with an 
anti-symmetric input profile, 

u{-l, y) = u(0, y) = 0 

u(x, 0) = — 2x — 1 on x = ( — 1,0) 

The exact solution to this problem contains symmet- 
ric expansion fans on the sides and a compression fan 
at the inflow that coalesces into a vertical shock at 

(•*5 y ) ( — 2 ’ 2 ) ■ 

The first mesh is cut-cartesian containing 26 x 26 
nodes. The FS and FV solutions, both using the Min- 
mod limiter, are presented in Figs. 14 and 15. Both 
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Contour 

spacing: 

0-1 

0.1 increments 


Fig. 14 Fluctuation splitting, Burgers equation. Fig. 16 Fluctuation splitting, Burgers equation, 

absolute error. 



Contour 

spacing: 

0-1 

0.1 increments 


Fig. 15 Finite volume, Burgers equation. 


Fig. 17 Finite volume, Burgers equation, absolute 
error. 


algorithms exhibit the same grid dependence on the 
amount of artificial dissipation as seen before, with 
the left-half solutions having more diffusion than the 
right halves, due to the grid orientation. Both meth- 
ods perform the same in the compression-fan region, 
coalescing into a shock to within the accuracy of the 
input-profile discretization. 

The shock is more sharply defined by FS than by 
FV. Fig. 14 has the correct shock speed, with nearly 
the entire gradient captured in one cell thickness. In 
contrast, Fig. 15 shows a slightly incorrect shock speed 
when using FV, as the shock progresses to the left 
beyond the coalescence point, even though the dis- 
cretization is conservative. The incorrect shock speed 
results from a non-symmetric distribution of the de- 
pendent variable to the left and right of the shock, 
caused by the excessive artificial diffusion generated 
on the grid-misaligned (left-hand) side. 

Contours of the absolute value of the error are pre- 
sented in Figs. 16 and 17. Errors from both computed 


solutions show a lack of symmetry, again reflecting the 
grid dependence of the artificial diffusion terms. The 
error levels from FS are less than from FV. The shock 
curvature in the FV solution at the coalescing point is 
clearly visible in Fig. 17, resulting in significant down- 
stream errors in the shock location as compared with 
the FS errors. 

This problem is repeated on a 25 x 25 mesh with 
symmetric diagonal cuts, favorably aligned with the 
advection directions. The FS and FV solutions, 
Figs. 18 and 19, are in good agreement. Plots of 
the absolute error contours, Figs. 20 and 21, show FS 
to be a little more accurate than FV for this case. 

The final mesh for this case is a truly unstructured 
triangulation containing 847 nodes and 1617 cells. The 
nodes are clustered to the outflow boundary, with a 
bias towards the left-hand side. The FS solution is 
presented in Fig. 22, showing very accurate and crisp 
shock resolution and good symmetry in the solution 
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Contour 

spacing: 

0-1 

0.1 increments 


Fig. 18 Fluctuation splitting, Burgers equation, Fig. 21 Finite volume, Burgers equation, absolute 
symmetric mesh. error. 



U = -2x -1 X 


Fig. 19 Finite volume, Burgers equation, symmet- Fig. 22 Fluctuation splitting, Burgers equation, 
ric mesh. unstructured mesh. 




Fig. 20 Fluctuation splitting, Burgers equation, Fig. 23 Finite volume, Burgers equation, unstruc- 
absolute error. tured mesh. 
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Fig. 24 Pure-diffusion problem error, diffusion Fig. 25 Pure-diffusion problem error, diffusion 
terms from Eqn. 8. terms from Eqn. 26. 


contours despite the mesh-clustering bias. In contrast, 
the FV solution in Fig. 23 has a more-diffuse shock and 
again an incorrect shock speed, with the outflow shock 
offset to the left of x = — A. The FV solution is also 
somewhat less symmetric than the FS solution. 

Linear Diffusion 

Choosing F = 0, the heat-conduction equation is 
obtained from Eqn. 1, 

Ut = V • (vWu) 

The test problem, a steady-state boundary value 
problem on a unit square, is taken from Tomaich. 4 
The Dirichlet boundary values are, 

u{-l,y) = 0, u(0,y) = sm(iry) 

u(x, 0) = 0, u(x, 1) = — sin(7r;c) 

The analytical solution on x = [—1,0], y = [0, 1] is, 

u(x,y ) = — t — [sinh(7r(;c + l)) sin(7ry) 
smh 7 r 

-F sinh(7ry) sin(7r(;r + l))] 

Both diffusion discretizations, Eqns. 8 and 26, are 
compared on a 438-node unstructured mesh. Fig- 
ures 24 and 25 plot the absolute value of the error 
in the converged solutions using Eqns. 8 and 26, re- 
spectively. 

The finite element treatment is clearly more ac- 
curate, and is used to discretize the diffusion terms 
for both FV and FS in the following section. The 
average-gradient results in Fig. 24 appear to exhibit a 
decoupling mode, similar to odd/even decoupling for 
structured meshes. 


Linear Advection-Diffusion 

The final test case is a linear advection-diffusion 
problem of Smith and Hutton. 18 The flux function 
is F = Xu, with, 

A = ( 22/(1 - .x 2 ), — 2x(l - y 2 )) 

The streamlines for this problem, while not truly circu- 
lar, are similar in orientation to the circular advection 
problem. The inflow profile is, 

u(x, 0) = 1 + tanh(20x + 10) 

The diffusion coefficient is chosen to be a constant, 
v = 10 -3 . No closed-form solution is known to this 
problem. 

A sequence of five unstructured meshes is consid- 
ered. The meshes, in order of increasing resolution, 
are: (A) 134 nodes, (B) 495 nodes, (C) 1928 nodes, 
(D) 7529 nodes, and (E) 28,915 nodes. The meshes 
have no preferred clustering or stretching and have 
nominal node-spacings of 0.1, 0.05, 0.025, 0.0125, and 
0.00625, labeled as Meshes A-E, respectively. 

Z/o-norms of the artificial and physical dissipations 
computed using both FS and FV are presented for each 
mesh in Table 1. Notice that the norm of the artifi- 
cial dissipation for FS is lower than the norm of the 
physical dissipation on Meshes D and E. In contrast, 
the artificial dissipation from FV is still larger than 
the physical dissipation even on the finest grid. Since 
the algorithms select only the larger of the physical or 
artificial dissipation (Eqn. 27), Table 1 suggests FS is 
grid resolved on Mesh D. 

Further evidence of a grid-resolved FS solution is 
seen in Figs. 26 and 27. The FS solution on Mesh E 
at the outflow boundary is presented along with the 
inflow profile and the corresponding pure-advection 
(v = 0) FS solution in Fig. 26. The pure-advection 
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Table 1 Lj-norms (xlO 6 ) of artificial and physical 
viscosities for advection/diffusion problem. 


mb 

(art.) 

FS 

\\4>v 2 

(phvs.) 

Mesh 

mb 

(art.) 

FV 

\\<f>v\\2 

(phvs.) 

1274 

215 

A 

1918 

190 

597 

265 

B 

640 

176 

192 

161 

C 

153 

107 

54 

76 

D 

123 

59 

13 

36 

E 

45 

34 



Fig. 26 Fluctuation splitting profiles on finest 
mesh, advection/diffusion problem. 


solution is seen to replicate the inflow profile ex- 
tremely well, with a clear separation from the diffused, 
v = 10 -3 , solution. Plotting only the FS results with 
respect to grid refinement, Fig. 27 shows a convergence 
of the outflow profile by Mesh C for FS. 

The accuracy of FS and FV are compared in Fig. 28, 
where the Mesh-E outflow solutions from FS and FV 
are plotted along with the Mesh C FS result. Taking 
the FS Mesh E solution to be the “truth” solution, it 
is clear that FS on Mesh C is more accurate in this 
case than FV is on the finest mesh. 

Computational efficiencies of the two algorithms are 
compared in Fig. 29, where the To-norm of the resid- 
ual is plotted versus CPU time for the fine-mesh FS 
and FV solutions, along with the FS convergence his- 
tory on Mesh D. The Mesh-E FS solution converges in 
700 sec on a 195 MHz SGI R10000 CPU. The corre- 
sponding FV solution takes 2.7 times longer than FS, 
due, in part, to the need to reconstruct gradient infor- 
mation at each node with FV for second-order spatial 
accuracy. However, considering the solution time to 
reach a given accuracy, it is more reasonable to com- 
pare the FS solution time on Mesh D to the finest-mesh 



Fig. 27 Fluctuation splitting grid convergence, ad- 
vection/diffusion problem. 



Fig. 28 Fluctuation splitting and finite volume for 
advection/diffusion problem. 

FV solution. The FS Mesh-D solution took only 80 sec, 
a factor of 23.6 times less than FV on Mesh E, and still 
shows better accuracy than the fine-mesh FV solution. 

An even greater speedup is seen with FS in con- 
junction with the van Albada limiter, where now the 
Mesh-B solution over-plots the curve from the finest 
grid, shown in Fig. 30. The corresponding FV result 
using the van Albada limiter on Mesh B is included, 
and clearly falls short of the FS accuracy. The FV 
case was repeated with the highly-compressive Super- 
bee limiter with little improvement in accuracy. The 
solution time for FS on Mesh B is three seconds, yield- 
ing a speedup factor of 600 over FV. 

The final set of results addresses convergence issues 
while pushing the positivity limits. Figure 31 com- 
pares two convergence histories for the second-order 
FS on Mesh B. The non-converging, though stable, 
convergence history is the result of using strict posi- 
tivity arguments to set the timestep (Eqn. 36). The 
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Fig. 29 Convergence histories for advec- 
t ion/ diffusion problem. 



Fig. 30 Advection/diffusion results using van Al- 
bada limiter. 

resulting solution is bounded and approximately cor- 
rect but oscillatory — attributed to limiter “ringing”. 
Full convergence is achieved by using first-order pos- 
itivity coefficients, which are not dependent on the 
limiters. The resulting local timesteps will not be as 
large as true second-order positivity would allow, but 
appear to be more robust. 

Summary of Results 

Fluctuation splitting and dimensionally-split finite 
volume schemes are compared in detail as applied 
to scalar advection, diffusion, and advection-diffusion 
problems. The fluctuation splitting scheme is seen to 
introduce less artificial dissipation while treating ad- 
vection terms, allowing for more accurate resolution of 
weakly dissipative advection-diffusion problems. The 
ability to resolve solutions to these problems on coarser 
meshes makes the fluctuation splitting scheme the pre- 
ferred choice over dimensionally-split finite volume. 



Fig. 31 Convergence rates using first- and second- 
order positivity coefficients. 

Linear advection test problems are utilized to inves- 
tigate the dependence of artificial diffusion production 
on grid orientation. Both fluctuation splitting and fi- 
nite volume are shown to exhibit grid dependencies, 
but with fluctuation splitting producing less artificial 
dissipation on all grids considered. 

A non-linear coalescing shock problem further ex- 
plores grid dependencies as cases are constructed that 
result in incorrect shock speeds for finite volume. Fluc- 
tuation splitting shows correct shock speeds for all 
grids and provides tighter shock capturing than finite 
volume. 

An advection-diffusion problem with small physical 
dissipation (diffusion coefficient of 10 -3 ) is considered 
where the reduction in artificial dissipation with fluc- 
tuation splitting results in a significant accuracy im- 
provement over finite volume. Convergence times are 
compared, showing a speedup of 2.7 for fluctuation 
splitting over finite volume on identical grids, using 
an explicit point Gauss-Seidel relaxation. However, 
a grid convergence study shows fluctuation splitting 
has better resolution of the solution on a coarser mesh 
than finite volume does on the finest mesh, resulting 
is a speedup of 23.6 for fluctuation splitting over finite 
volume. 

Based upon these significantly reduced solution 
times for solving model problems, as compared to 
the current state-of-the-art dimensionally-split finite 
volume method, fluctuation splitting is considered a 
worthwile scheme to pursue for modeling fluid dynamic 
problems. 
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